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Part  I  -  Technical  Manual 


I,)  Introduction 

The  development  of  all-eleotrie  aireraft  is  a  high  priority  in  the  avionies  eommunity.  Current 
aireraft  use  a  eombination  of  hydraulie,  pneumatic,  and  electric  systems  for  flight  control. 
However,  the  expectation  for  future  airplanes  is  a  single,  electric  system  using  electromechanical 
actuators  (EMAs).  Such  a  system  would  reduce  the  cost  to  build,  operate,  and  maintain  aircraft. 
It  would  also  make  aircraft  lighter,  more  reliable,  safer,  and  more  easily  reconfigurable,  reducing 
the  turnaround  for  new  technology. 

One  of  the  greatest  hurdles  to  replacing  all  hydraulic  actuators  with  EMAs  is  heat  generation,  a 
consequence  of  the  absence  of  cooling  hydraulic  fluid.  Accurately  quantifying  the  heat 
generated  is  complicated  by  the  highly  transient  and  localized  nature  of  the  power  demands  of  an 
EMA's  motor,  an  especially  significant  issue  in  aircraft.  Thus,  accurate  modeling  must  be 
dynamic. 

The  parameters  that  characterize  the  motor  are  nonlinear  functions.  In  steady-state  analysis, 
these  nonlinearities  can  simply  be  averaged  out,  and  often  the  parameters  are  considered  to  be 
constant  [1].  But  during  transient  analysis  these  nonlinearities  can  become  critical  [2]. 
Consequently,  time-averaged,  linear  models  [3]  are  inadequate  for  thermal  analysis  and 
management  designs.  So,  proper  treatment  of  the  problem  requires  dynamic,  nonlinear  analysis. 

The  motor  parameters  can  be  determined  by  careful  experimentation  [2-3],  but  this  is  a  time 
intensive  process  and  can  be  expensive  especially  when  considering  multiple  potential  designs 
for  a  motor.  Jens  Otto  suggested  using  a  reduced  order  coenergy  model  [4],  which  yielded 
results  comparable  to  EEM  results.  Eikewise,  Roshen  considered  more  involved  empirical 
formulas  which  accounted  for  excess  eddy  current  losses  [5]. 

Through  EEM  modeling,  the  nonlinear  parameters,  such  as  self  and  mutual  inductances  of  the 
three  phase  windings,  can  be  determined  [6-7].  Though  dynamic  EEM  can  be  very  accurate  [4], 
such  detailed  simulation  is  very  slow.  Eortunately,  EEM  is  only  necessary  to  quantify  the 
parameters;  it  is  not  necessary  for  simulating  the  motor  over  lengthy  mission  profiles.  Instead,  a 
nonlinear,  lumped  element  model  (NE-LEM)  can  use  the  parameters  from  the  EEM  model  and 
then  dynamically  simulate  the  motor  just  as  accurately  as  the  EEM  model  but  at  a  much  lower 
computational  cost.  To  our  knowledge,  we  are  the  first  to  address  the  nonlinear  dynamic 
modeling  of  a  permanent  magnet  motor  and  describe  both  the  control  and  thermal  performance 
of  the  motor  in  following  highly  transient  mission  profiles. 
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Figure  1.1  is  the  3-D  model  of  a  typieal  Permanent  Magnet  Synehronous  Maehine  (PMSM) 
servo  motor.  This  design  features  a  12-slot  stator  and  a  10-pole  rotor  (Figure  1.2).  In  the 
following,  we  will  diseuss  the  proeedures  for  eleetrieal  and  thermal  modeling. 


Figure  1.1.  3-D  model  of  the  PMSM  motor  design. 


Figure  1.2.  Cross-section  of  motor  and  rotor. 
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II,)  Simulation  Flow  Chart 

Before  describing  the  details  of  the  simulation,  we  give  the  overall  flow  chart  in  Figure 
1.3  as  follows; 


Figure  1.3;  Flow  chart  for  simulation 
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For  the  motor  electrical  simulation  block,  the  flow  chart  is  given  in  Figure  1.4. 


Figure  1.4;  Flow  chart  for  motor  electrical  simulation 


III.)  Electrical  Model 

1.  Dynamical  Equations 

The  motor  dynamics  are  modeled  by  four  primary  dynamical  equations  in  direct- 
quadrature  reference  frame  (dqO); 
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where  Uj  is  direct  input  voltage,  is  quadrature  input  voltage,  is  direct  current, 
is  quadrature  current,  is  phase  resistance,  is  direct  inductance,  is 
quadrature  inductance,  p  is  the  number  of  poles,  is  the  mechanical  frequency 

multiplied  by  the  number  of  pole  pairs  pl2  ,  is  the  flux  linkage  from  the 
permanent  magnet,  is  the  motor  torque  generated  by  the  magnetic  fields,  Tp  is  the 
load  torque,  is  the  mechanical  angular  acceleration  multiplied  by  the  number  of 
pole  pairs  p  12  ,  I  is  the  rotor's  moment  of  inertia,  and  c  is  the  rotor's  coefficient  of 
friction  from  windage  and  bearings. 

where  u  is  input  voltage,  is  direct  current,  is  quadrature  current,  is  phase 
resistance,  is  direct  inductance,  is  quadrature  inductance,  p  is  the  number  of 
poles,  is  the  mechanical  frequency  multiplied  by  the  number  of  pole  pairs  pl2  , 
Apj^  is  the  flux  linkage  from  the  permanent  magnet,  is  the  motor  torque  generated 
by  the  magnetic  fields,  Tp  is  the  load  torque,  is  the  mechanical  angular 
acceleration  multiplied  by  the  number  of  pole  pairs  p  12  ,  I  is  the  rotor's  moment  of 
inertia,  and  c  is  the  rotor's  coefficient  of  friction  from  windage  and  bearings. 

The  motor  losses  are  related  to  the  motor  parameters,  such  as,  R^,  c ,  L^,  and  . 

Since  losses  have  an  effect  on  motor  behavior,  they  should  be  modeled  through  the 
motor  parameters,  not  calculated  in  post-processing,  as  is  often  done  [5].  The  copper 
loss  can  be  calculated  as 

Pcu  =  ka  +  ’  (3) 

or  equivalently 

(6) 
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for  balanced  3-phases.  In  either  case,  note  that  the  resistance  is  the  parameter  directly 
associated  with  the  power  loss  in  the  windings.  Likewise,  windage  and  bearing  losses 
are  directly  associated  with  the  coefficient  of  friction,  c  .  Ail  iron  losses  (hysteresis, 
classical  eddy,  and  excess  eddy  losses)  can  be  associated  with  the  nonlinear 
inductance,  L . 

Much  of  the  motor  modeling  researeh  done  in  the  past  [1]  and  [8]  and  even  recent 
work  [9]  use  constant  parameter  values  in  simulations,  although  some  did  use 
nonlinear  inductances  [2]  and  [7],  and  [10].  Since  the  BH  curve  is  nonlinear  and  H 
is  a  function  of  current,  inductance  is  also  a  nonlinear  function  of  current.  For  our 
model,  we  used  FEM  to  obtain  and  as  functions  of  .  At  this  point  we  are 

neglecting  changes  in  since  our  control  algorithm  already  seeks  to  maintain  a  zero 
value  for  .  The  inductance  curves  we  do  have  were  calculated  in  ANSYS  using  the 

magnetization  curve  for  our  particular  soft-magnetic  material  and  the  geometry  of  our 
motor.  The  analysis  was  done  for  zero  direct  current  and  for  a  wide  range  of 
quadrature  currents  (Figure  1.5). 


Inductance  as  a  Function  of  Current 


Figure.  1.5.  Direct  and  quadrature  inductance  as  a  function  of  current. 


Time  Stepping 

The  size  of  variable  arrays  affects  the  speed  of  a  simulation.  Using  prerecorded  input 
profiles,  high  resolution  is  necessary  to  accurately  capture  transient  moments.  A 
uniform,  high  resolution  would  result  in  a  very  large  profile.  To  keep  our  array  sizes 
small  but  maintain  high  resolution  where  needed,  we  used  a  line  simplification 
algorithm  to  reduce  the  number  of  points  where  little  activity  would  occur.  This 
compression  technique  allowed  us  to  reduce  the  profiles  to  less  than  4%  of  their 
original  size. 
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That  said,  the  electrical  time  constant  is  often  much  smaller  than  the  time-step  size  of 
the  compressed  profiles.  Since  the  frequency  of  the  electrical  component  of  the 
simulation  is  directly  related  to  the  first  derivative  of  the  input  angular  position 
profile  for  the  motor,  we  can  dynamically  scale  the  electrical  time  step  non- 
uniformly.  Therefore,  we  had  two  loops,  one  macro  loop  that  stepped  through  the 
elements  of  the  profile  and  a  nested  micro  loop  with  very  small,  but  variable  time 
steps  for  the  electrical  equations.  The  result  of  this  apparent  complication  was  a 
many-fold  speed  improvement.  This  is  a  tremendous  advantage  of  not  relying  on 
MATLAB  built-in  ODE  solvers,  which  do  not  have  the  advantage  of  knowing  how 
big  of  a  time  step  to  take  but  must  rather  continually  search  for  it. 


Field  Oriented  Control  (FOC) 

A  thorough  analysis  of  motor  heat  generation  should  include  how  the  motor  is  driven. 
Our  control  equations  are  derived  from  the  primary  simulation  equations  (1)  through 
(4).  Starting  with  equation  (4)  and  noting  that  in  (3)  can  be  factored  out,  we  get 
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which  can  be  rearranged  to  be 


h  = 
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At  this  point  we  reinterpret  which  is  dco^^  /  dt  to  be  I  dt ,  where  is  the 
error  in  ,  with  a  tuning  coefficient,  ,  in  front.  This  gives 
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where  is  the  desired  current  for  the  next  time  step  and  is  the  desired  next- 
step  angular  speed  equal  to  .  Using  (4)  to  substitute  for  in  (9),  we 

arrive  at  (10): 
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Recalling  that  {co^^  -  equals  ,  (10)  simplifies  to 
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This  will  track  the  desired  motor  velocity  well,  but  there  will  be  some  rotor  angle  drift 
over  time.  To  counter  this  drift  error,  we  added  a  theta  error  term,  k^eg  /  dt : 


i„  = 
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dt 
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(12) 


The  kp  and  k.  coefficients  are  the  only  PI  values  that  need  to  be  tuned  in  the  control 

algorithm.  For  maximum  efficiency,  it  is  generally  desired  that  the  direct  current  be 
zero  since  any  flux  generated  in  the  direct  axis  would  not  contribute  to  torque. 
However,  direct  current  does  not  necessarily  directly  translate  to  direct  flux.  As  the 
following  equation  shows,  direct  current  can  contribute  to  useful  torque  if  the  direct 
and  quadrature  inductances  are  not  equal,  as  given  in  (3). 

Most  of  the  time,  however,  the  desired  direct  current  set  point,  i* ,  should  be  zero. 
With  the  desired  currents  determined,  the  desired  voltages  can  be  found; 
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(14) 


u  =  R  i  +  L 

q  s  q  q 
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dt 


^PM  ’ 


^  * 

where  and  are  based  on  the  desired  eurrent  values.  It  should  be  noted  that 

while  many  of  the  variables  are  updated  based  upon  the  desired  currents,  was  not 
updated  for  the  control  equations.  This  is  a  reasonable  simplification  since  for  a 
round  rotor  is  independent  of  should  normally  be  zero,  and  the  factor 

is  divided  out: 


—  =  (15) 

For  the  non-round  rotor  case,  our  direct  interpolation  method  could  be  used  to  include 
the  dependence  on  theta. 

Even  after  the  desired  control  values  were  calculated,  the  reality  of  limits  (current  and 
voltage)  needed  to  be  included.  Instead  of  directly  capping  those  values,  we  realized 
that  the  system  converged  more  quickly  and  ran  more  smoothly  if  the  functions  used 
were  all  smooth  [11].  For  our  case,  we  used  a  segmented  curve  of  first  and  second- 
order  polynomials. 


IV.)  Thermal  Model 
1,  Model  Setup 

A  lumped  node  thermal  network  is  to  represent  the  temperature  of  every  solid  part  of 
the  EMA.  The  thermal  resistances  and  capacitances  between  the  nodes  can  be  treated 
as  electrical  resistances  and  capacitances.  Hence  the  temperature  of  every  node  can  be 
solved  as  the  voltage  in  this  equivalent  network  (Table  1.1).  This  approach  is  well 
developed  in  motor  design  industry  [12,  13].  Some  commercial  motor  design 
software  package  has  already  included  the  thermal  network  simulation,  i.e.  Motor- 
CAD  [14].  Some  studies  have  been  made  with  FEA  analysis  and  experimental  testing 
have  shown  that  such  an  approach  is  valid  [15]. 

Table  1.1  The  analogy  of  equivalent  thermal  circuit 
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Equations  (16)  and  (17)  show  how  to  calculate  of  the  R  and  C  values  for  simple  one¬ 
dimensional  geometry.  Aeeurate  estimates  for  R  and  C  values  are  not  possible  for 
eomplieated  geometries  sueh  as  an  eleetrie  motor.  The  traditional  lumped  node 
thermal  network  is  based  on  a  “forward”  direetion  modeling  proeess.  In  this  process, 
with  all  the  material,  geometry  and  eonstruction  information  provided,  the  thermal 
resistanee  R0  and  eapacitance  C0  are  ealculated  based  on  empirical  and  theoretical 
relations.  For  example,  when  the  winding  wire  type,  winding  structure,  potting 
material,  slot  width  and  teeth  thickness  are  given,  the  thermal  resistance  from  winding 
to  stator  and  from  winding  to  rotor  can  be  calculated.  Because  this  type  of  approach 
requires  large  number  of  empirical  relations  to  aehieve  high  accuracy,  a  “reverse” 
modeling  process  is  introduced  in  this  paper.  A  detailed  3-D  solid  model  of  a  target 
motor  is  constructed  and  imported  to  an  FEA  software  like  ANSYS  [16]  to  perform 
steady-state  and  transient  simulation.  With  these  results  we  can  estimate  and  select 
the  values  for  the  thermal  resistances  and  capacitances.  The  advantage  of  this  method 
is  that  we  can  evaluate  the  fidelity  of  the  lumped  node  model  with  a  real  motor,  add 
or  reduce  nodes  to  inerease  the  model  aceuraey  and  effieieney,  and  estimate  the 
maximum  error  between  the  node  temperature  and  maximum  temperature  in  real 
motor. 


This  procedure  ean  also  be  extended  to  model  the  gear-box,  motor-driver  and  drive- 
train,  and  even  include  the  aircraft  wing  surfaces  and  frames.  The  thermal  network 
ean  also  be  ineorporated  into  a  multi-physios  model  to  simulate  the  eleotrioal,  thermal 
and  meohanioal  performance  of  the  whole  EMA  and  its  supporting  struoture.  Once 
the  proper  thermal  resistanoes  and  oapacitanoes  are  seleoted,  this  simulation  engine 
can  be  used  with  various  time  dependent  boundary  oonditions  during  the  whole 
mission  duration,  inoluding  the  air  temperature,  airoraft  speed,  altitude  and  sunlight. 
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The  computational  requirement  of  such  a  simulation  is  negligible  compared  to  the 
FEA  simulation. 

Due  to  the  symmetry,  the  geometry  of  the  model  shown  in  Figure  1.1  is  divided  into  a 
quarter  section  and  imported  into  ANSYS  for  the  FEA  simulation.  The  nodes  are 
numbered  in  Figure  1.6.  These  numbers  are  also  the  same  as  those  in  the  lumped 
node  network  shown  in  Figures  1.7. 


Air  gap  (Node  10) 

Magnet  (Node  9) 

Rotor  core 
(Node  6) 

Case 
(Node  4) 

Case  surface 
(Node  5) 

Stator  iron 
(Node  3) 

Copper  winding  (Node  1) 


Figure  1.6.  Thermal  node  locations  on  motor. 
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Detailed  Lumped  node  model  Of  EMA  (transient) 


I.  copper 

3.  stator 

4.  case  and  back  cover 

5.  case  surface 

6.  shaft  body 

7.  shaft  front  end 

8.  shaft  back  end 

9.  magnet 

10.  air  gap 

II.  connection  gear 

12.  front  bearing 

13.  back  bearing 

14.  front  cover 

15.  moving  mechanism 

16.  wing  surface 

17.  gear  box  case 


R16 


Airflow  outside 


Figure  1.7.  Lumped  node  model  of  EMA 

2,  Dynamic  Thermal  Equations 

After  all  the  thermal  resistor  and  capaeitor  values  of  the  lumped-node  network  model 
of  the  motor  in  Figure  5  are  known,  the  same  values  of  the  thermal  resistors  and 
eapaeitors  ean  be  used  to  simulate  the  temperature  response  of  the  motor  parts  with 
any  combination  of  heat  losses  and  boundary  conditions.  For  this  standard  resistor- 
capacitor  electrical  network,  a  set  of  first-order  ordinary  differential  equations  can  be 
written  as; 


U.-U,  ,  ^10- 

7?3  R 

(7?g  )f/jQ 


Uy-U,  ,  ^  dU^ 

R,  ’  dt 


f/3  f/i-f/3  ^  dU^ 
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dt 


-U^Rg  +RloUg  +RgR^Qli^{t) 


(18) 

!A>) 

(19) 

(20) 

12 


U,-U,  dU,  _ 

Ret  R,  "  dt 

i?13  R^  dt 


(21) 

(22) 


^  ^  ^,-^3  U,-U, 

R^a  R\7,  R^  Ra 

Un  -U,,  ,  U,  -U,,  ^  dU,,  _ 

I  ^12  “ 

Ri2  Rj  dt 

U,-U,  ^10-^9,^  dU,  _ 

i?g  Rg  ^  dt 

Ue-Ui  Ue-U,  ^9-^6,^ 
^6.  Ret  Re 

U^-Un  ,r  ^  = 

R,  Re.  '  dt 

Un-Uu  Uu-Ua.^  dU,, 
Rn  Ra.  dt 


^  BF  (0 


^  M  (0 


dt 


-I 


Rla 


=  1 


7?14 


=  0 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


Because  the  time  constant  (about  10  ms)  for  the  thermal  component  is  much  larger 
than  for  the  electrical  component  in  our  case,  the  thermal  component  was  calculated 
less  often,  thereby  increasing  the  speed  of  the  simulation. 
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Part  II  -  User's  Manual 


I,)  Introduction 


a.  Purpose  of  the  Software 

This  software  provides  integrated  non-linear  dynamic  modeling  including  both 
electrical  model  and  thermal  model  together  with  a  field  oriented  control  scheme. 
Different  motor  configurations  can  be  entered  into  the  software.  The  simulation 
provides  the  user  with  all  the  necessary  motor  information,  such  as;  mission  profile 
following,  motor  torque,  motor  current,  phase  voltage,  input  power,  power  losses  in 
the  windings,  and  a  thermal  profile.  This  comprehensive  and  detailed  look  at  motor 
control  together  with  heat  generation  and  transfer  will  provide  solid  foundation  for 
users  to  design  highly  efficient  EMAs. 

b.  Software  Capability 

Package  of  a  nonlinear  dynamic  modeling  of  a  permanent  magnet  motor  that 
describes  both  the  control  and  thermal  performance  of  the  motor  in  following  highly 
transient  mission  profiles.  One  of  the  most  attractive  features  of  this  model  is  that  it 
is  able  to  incorporate  various  motor  designs.  The  software  manual  assumes  that  the 
mission  profile  is  provided  to  the  user.  For  users  familiar  with  control  coding  will  be 
able  to  manipulate  some  constants  in  the  control  helping  to  tune  for  their  specific 
motor  designs. 


c.User  Privileges 

The  main  advantages  to  the  user  is  to  put  in  their  own  motor  parameters  for  alternate 
designs  and  test  whether  the  design  is  capable  of  providing  proper  control  and  also 
find  out  the  thermal  performance.  The  motor  parameter  code  (located  in  Appendix 
A)  contains  a  variety  of  values  that  must  be  determined  by  the  user  if  an  alternate 
design  is  used.  Values  such  as  phase  resistance,  number  of  magnetic  poles,  and  the 
moment  of  inertia  must  be  determined.  Motor  limitations  are  also  set  in  this  file  to 
ensure  that  these  limits  are  not  exceeded  during  operation.  The  inductance  values,  as 
well  as  the  thermal  resistance  and  capacitance  values,  can  be  obtained  from  FEM 
models.  If  a  new  motor  model  is  to  be  used  that  maintains  the  non-linear  nature  of  the 
simulation,  the  new  inductance  values,  thermal  resistance  and  capacitance,  should  be 
recalculated  using  FEM.  Although  determining  these  values  is  outside  the  scope  of 
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this  manual,  the  user  is  able  to  manipulate  any  of  the  data  inside  the  motor  simulation 
to  suit  the  needs  of  their  design. 


II,)  Input  Parameters 


There  are  the  two  file  names  ‘motorParameter’  and  ‘missionProfile.xlsx’  for  input 
motor  parameters  and  mission  profiles. 

motorParameters.m  eontains  motor  parameters  speeifie  to  a  partieular  motor  design  as 
well  as  ealeulated  parameters.  If  you  are  to  ehange  the  name  of  this  folder  or  make  a 
new  motor  design,  the  new  name  of  the  folder  must  be  input  into  the  main  program. 
The  code  is  displayed  in  Part  IV  but  these  parameters  show  the  motor  parameters  for 
this  particular  example.  These  values  can  be  changed  to  accommodate  other 
motor/EMA  designs. 

Tables  2. 1-2.4  summarize  the  input  parameters  in  motorParameters.m. 


Table  2.1.  Motor  Electrical  Parameters 


Motor 

Electrical 

Parameters 

Explanations 

Units 

Source 

Rs  room 

phase  resistance  at  room  temperature 

Q 

EEM  or 
Measurement 

P 

number  of  poles 

Given 

I 

inertia  moment  of  rotor 

kgm^ 

Given 

c 

friction  coefficient 

kg  m  /s 

Given 

lambdaPM 

permanent  magnet  flux  amplitude  kpM 

Wb 

EEM  or 
Measurement 

EdSN 

array  of  direct  inductance  Ed 
corresponding  to  current  array  iq 

H 

EEM  or 
Measurement 

EqSN 

array  of  quadrature  inductance  Eq 
corresponding  to  current  array  iq 

H 

EEM  or 
Measurement 

iq 

array  of  quadrature  current 

A 

Given 
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Table  2.2.  Motor  Thermal  Parameters 


Motor  Thermal 
Parameters 

Explanations 

Units 

Source 

Rth 

thermal  resistance 

T/W 

FEM 

Cth 

thermal  capacitance 

J/°C 

FEM 

Troom 

room  temperature 

°C 

Given 

IC 

winding  loss 

w 

Electrical 

Simulation 

IS 

stator  loss 

w 

Electrical 

Simulation 

IW 

winding  loss 

w 

Electrical 

Simulation 

IBB 

rear  bearing  loss 

w 

Mechanical 

Simulation 

IBF 

front  bearing  loss 

w 

Mechanical 

Simulation 

IM 

magnet  loss 

w 

Electrical 

Simulation 

IR7a 

conduction  heat  loss  to  gear  box  axle 

w 

FEM 

IRI4 

conduction  heat  loss  to  gear  box  case 

w 

FEM 

Table  2.3.  Drive  Train  Parameters 


Drive  Train 
Parameters 

Explanations 

Units 

Source 

Ncr 

total  coupling  ratio 

rad/m 

Given 

eflDrive 

drive  train  efficiency 

Given 

theta  meintial 

initial  rotor  angle  multiply  p/2 
when  stroke  is  0 

rad 

Given 

Table  2.4.  Simulation  Limits  (These  parameters  are  given.) 


Simulation  Eimits 

Explanations 

Units 

xMinEim 

minimum  value  of  stroke 

m 

xMaxEim 

maximum  value  of  stroke 

m 

tauMEim 

motor  torque  limit 

Nm 

uEim 

voltage  limit 

V 

iEim 

current  limit 

A 

didtEim 

current  change  rate  limit 

A/s 

vEim 

speed  limit 

m/s 
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Excel  file  missionProfile.xlsx  provides  data  for  mission  profile.  The  mission  profile 
data  is  assumed  to  have  been  provided  already  and  is  not  subjeet  for  diseussion  in  this 
manual.  The  input  parameters  in  mission  profile  is  summarized  in  Table  5. 


Table  2.5.  Mission  Profile  Parameters 


(may  eome  from  aerodynamie  measurement  or  simulation) 


Mission  Profile 
Parameters 

Explanations 

Units 

tN 

array  for  time 

s 

xSN 

array  for  stroke  profile 

M 

FEN 

array  for  load  foree  profile 

N 

Due  to  the  faet  that  there  can  be  various  mission  profiles  and  motor  designs  it  is 
important  to  ensure  that  your  file  names  in  the  ‘FOC’  folder  mateh  the  file  names  as 
represented  in  the  eode  above. 

There  are  also  input  parameters  used  to  set  up  simulation.  These  parameters  are 
summarized  in  Table  2.6. 


Table  2.6.  Simulation  Parameters  (These  parameters  are  given.) 


Simulation  Parameters 

Explanations 

Units 

spe 

samples  per  eyele  of  highest  frequeney 

tauTh 

thermal  time  constant 

s 

kp 

proportional  parameter  for  PI  controller 

ki 

integral  parameter  for  PI  controller 

A  s/rad 

alphadt 

forward  time-stepping  weight  for 
implieit  algorithm  in  eleotrieal  loop 

epsilon 

error  for  eleotrieal  oonvergenoe  loop 

Ncvg 

maximum  number  of  convergence 
iterations 

If  the  experience  of  the  user  is  adequate  in  the  design  of  the  control  scheme,  there  are 
few  parameters  that  ean  be  changed  to  deerease  error  in  the  ability  for  the  motor  to 
follow  the  flight  profile  data.  Users  ean  ehange  the  proportional  eonstant,  or  kp  value, 
and  the  integral  eonstant,  or  ki  value  of  the  eontrol. 
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%  Set  PI  coefficients. 


38 

39 


%  [] 

;  %  [A  s/rad] 


221 

222 

223  - 

224  - 

225 

226 

227 

228  - 

229 

230 

231  - 

232  - 

233  - 

234 

235 

236  - 

237 

238 

239 

240  - 


%  Get  iqS. 

tauM_iq  =  3*p/4* (lambdaPM  +  id* (Ld-Lq) ) ;  %  [N  m/A] 
iqS  =  iq  +  2/  (p*tauM_iq)  *(I/()cp*e’toega_me/dtE  -  alpha_me) 
(icp^Omega  me)  )  me/dtE;  %  [A] 


%  Set  ids . 

ids  =0;  %  id  should  almost  always  be  zero  [A] 


%  Rate  limit  iqS  (reduces  transients) . 

diqdtS  =  (iqS-iq) /dtE;  %  Predicted  current  rate  [A/s] 

diqdtS  =  sicrn (diqdtS) *min ( [didtLim, abs (diqdtS) ] ) ;  %  Cap  diqdt  [A/s] 

iqS  =  iq  +  diqdtS*dtE;  %  Rebuilt  iqS  [A] 


%  Saturation  limit  abs  of  iqS. 

iqS  =  sign (iqS) *( (abs (iqS) <i3) * (abs (iqS)  -  (abs (iqS) >=il) * . . . 
( (abs (iqS) -iLim* (1-rSat) ) "2) / (4*iLim*rSat) )  +  ... 

(abs (iqS) >=i3) *iLim)  ;  %  [A] 


:ga  meP  =  omega 


33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 


Two  parameters,  alpha  dt  and  beta  dt  located  in  the  main  code  are  weighted  averages 
used  in  calculating  next  step  values.  Generally  the  beta  value  is  set  to  .5  but  can  be 
altered  at  the  user’s  discretion  for  algorithm  convergence.  Another  parameter  that 
can  be  changed  to  help  improve  design  is  spc  or  samples  per  frequency  as  shown 
below. 

% - 

%  Set  up  simulation  parameters. 

% - 

^^^pc"  =  10;  N.  %  Samples  per  cycle  of  highest  frequency 

V^a^h  =  %  Thermal  time  constant  [s] 

%  Set  PI  coefficients. 

]tp  =  0.1;  %  [] 

kl  =  0.0001;  %  [A  s/rad] 

%  Set  differential  time-stepping  weights. 
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III.)  Software  Startup 


1 .  Once  all  the  proper  m-files  are  verified  click  on  the  ‘Start  Menu’  and  open  the 
MATLAB  software  through  the  program  listing. 

2.  Once  MATLAB  is  open  go  to  ‘File’  and  click  on  ‘Open’ 

3.  Go  to  through  the  above  mentioned  directory  path  to  find  the  ‘FOC’  folder 

4.  Click  on  the  ‘FOC’  folder  and  then  click  on  the  file  name  ‘FOC.m’ 

5.  The  following  code  should  be  displayed  in  the  MATLAB  editor  screen 


B  Editor  -  C\Users\WUTX\Desktop\Wor1dng\focm 


File  Edit  Text  Go  Cell  Tools  Debug  Desktop  Window  Help 


^  iff ' 

M  ^  1 

Stack:; 

Q  t3 


1.0 

+ 

•r 

1.1 

X 

%" 


1 

2 

3 

4 

5 
€ 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


%  Nonlinear  Dynamic  Modeling  and  Field  Oriented  Control 
%  of  Permanent  Magnet  (PM)  Motor 
% 

%  Written  in  SI  or  MKS  Dnit  System.. 

% 

%  Authors : 

%  David  Woodburn 

%  Dr.  Lei  Zhou 

%  Dr .  Thomas  X .  Wu 

clear  all;  close  all;  clc; 


%  Input  m.otor  param.eters. 

% - 

motor  Par  am.eter 
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IV,)  Software  Execution  and  Procedure 


1 .  In  order  to  run  the  simulation  you  can  either  click  in  the  editor  window  and  press 
the  F5  shortcut  key,  or  click  the  run  button  at  the  top  of  the  editor  window. 


Editor  -  0\Users\WUTXtDesktop\Wof1dng\focjn 


Frie  Edit  Terf  Go  Cell  Tools  Debug  Desktop  Window  Help 


D  ^ 


<)6 


•0  O'  I  ’  #4  ^  ^ 


«  «  V  a 


Stc 


-  1.0 


+  1.1 


1  %  Nonlinear  Dynaitic  Modeling  and  Field  Oriented  Control 

2  %  of  Permanent  Magnet  (PM)  Motor 

3  % 

4  %  Written  in  SI  or  MKS  Unit  System. 

5  % 


2.  The  following  menu  will  pop  up  to  add  the  foc.m  file  to  the  MATLAB  path. 


MATLAB  Editor 


J  j  File  C:\Users\WUTX\Desktop\Working\foc.m  is  not  found  in  the 
'  r  current  directory  or  on  the  MATLAB  path. 

To  run  this  file,  you  can  either  change  the  MATLAB  current  directory  or  add  its 
directory  to  the  MATLAB  path. 


Change  Directory 

Add  to  Path 

Cancel 

Help 

•  Click  Change  Directory.  Different  versions  of  MATLAB  will  show  different 
popup  boxes  when  the  directory  is  not  specified.  Simply  ensure  that  all  the  files 
needed  to  run  the  program  are  located  in  the  same  folder  before  changing  the 
path. 

3.  The  simulation  will  begin  running.  The  runtime  will  vary  based  on  the  system. 
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There  will  be  a  statistical  report  and  a  total  of  9  graphs  displayed  as  well  once  the 
simulation  is  finished. 


V,)  Output  Parameters 


The  Output  parameters  are  summarized  in  Tables  2.7  and  2.8  for  reporting  statistics  and 
plot  figures. 


Table  2.7.  Output  Parameters  to  Report  Statistics 


Parameters 

Explanations 

Units 

PcuMean 

mean  winding  copper  loss 

W 

tauMMax 

maximum  motor  torque 

Nm 

FMax 

maximum  actuator  force 

N 

vMax 

maximum  actuator  velocity 

m/s 

aMax 

maximum  acceleration 

m/s^ 

dtMax 

maximum  macro  time  step 

s 

dtEMax 

maximum  electrical  time  step 

s 

dtEMin 

minimum  electrical  time  step 

s 

Table  2.8.  Output  Parameters  in  Figures 


Output  Parameters 

Explanations 

Units 

tNR 

array  for  output  time  sequence 

s 

xN 

array  for  actual  stroke 

m 

tauMN 

array  for  motor  torque 

Nm 

udN 

array  for  direct  voltage 

V 

uqN 

array  for  quadrature  voltage 

V 

idN 

array  for  direct  current 

A 

iqN 

array  for  quadrature  current 

A 

PcuN 

array  for  instantaneous  copper  loss 

W 

PinN 

array  for  instantaneous  input  power 

W 

TthN 

2D  array  for  temperature  at  different 
nodes 
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VI,)  Expected  Results 


The  simulation  has  fully  run  when  the  following  statistics  are  displayed  in  the  command 
window  to  include; 

Mean  winding  power  loss,  max  motor  torque,  max  actuator  force,  max  velocity,  max 
acceleration,  min  electrical  time  step,  max  electrical  time  step  and  the  total  elapsed  time 
for  the  program  to  run. 


Command  Window 


Dynamic  Modeling  and  Field  Oriented  Control 
Mean  winding  power  loss  =  1.8149  W 
Max  motor  torque  =  2.1564  N  m 
Max  actuator  force  =  2873.9049  N 
Max  velocity  =  0.092486  m/ s 
Max  acceleration  =  2.685  m/s''2 
Min  electrical  time  step  =  0.0005  s 
Max  electrical  time  step  =  0.0039899  s 
Elapsed  time  is  300.858010  seconds. 


Also  figures  are  plotted  are  a  result  of  the  flight  mission  profile  used  specific  to  this 
manual  and  will  vary  based  on  other  motor  parameters  and  flight  mission  data.  This 
profile  was  a  five  minute  section  of  a  full  flight  profile.  Acceptable  values  will  be 
specific  each  set  of  parameters  or  profile  data;  however,  some  generalizations  can  be 
made  about  parameters  such  as  power  loss. 
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1.  The  stroke  profile  shows  that  movement  of  the  shaft  in  the  EMA.  The  green  line 
represents  the  aetual  movement  of  the  aetuator  based  on  the  mission  profile.  The  blue 
line  is  the  desired  movement  of  the  EMA.  ft  is  hard  to  see  the  blue  beeause  the  green 
is  eovering  indieating  good  eontrol.  The  desired  stroke  is  in  exeel  file 
missionProfile.xlsx.  The  aetual  stroke  is  saved  in  strokeAetual.xlsx. 


Mechanical  Stroke  Following 


Eigure  2.1  (a);  EMA  stroke 


To  further  see  the  aeeuraey  of  the  stroke  following  plot,  a  zoomed  seetion  of  Eigure  1 
is  shown  below. 


24 


Mechanical  Stroke  Following 


Figure  2.1  (b);  Zoomed  section  of  stroke 

2.  In  the  following  Figure  2.2  (a)  and  (b),  the  load  force  from  mission  profile  and  load 
torque  converted  by 


F, 


= 


Load 


^crldriv, 


where  Fcr  is  coupling  ratio  and  ridrive  is  drive  train  efficiency.  In  Figure  2.2  (c),  the 
magnetic  torque  is  shown.  Comparing  Figure  2.2  (b)  and  (c),  we  found  they  are  close. 
This  is  because  the  moment  of  inertia  of  the  motor  we  used  is  small.  The  magnetic 
torque  and  load  torque  are  related  by; 


J 


d 

dt^ 


The  magnetic  torque  is  saved  in  Excel  file  tauM.xlsx. 
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Lead  Torque,  r  [Nm] 


Load  Force 


Figure  2.2(a):  Load  Force 

Load  Torque 


Figure  2.2(b):  Load  Torque 
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Magnetic  TotqLe,  z:  [Nm] 


2.5 


Magnetic  Torque 


2 

1.5  - 
1  - 
0.5 
0  - 
-0.5  - 
-1 

-1.5- 

-2 

_2  5 _ ^ ^ ^ ^ ^ ^ ^ ^ ^ _ 

■  0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 

Time,  f  [min] 


Figure  2.2(c):  Magnetic  Force 


3.  The  dq  currents  provided  to  the  motor  are  represented  in  Figure  2.3.  Notice  that  the 
current  id  nearly  constant  around  0  and  doesn’t  exceed  the  dotted  blue  line  limitations. 
This  is  a  proper  reading  as  direct  current  does  not  contribute  to  the  rotation  of  the 
rotor.  Only  iq  contributes  to  the  torque  as  seen  by  the  variation  in  the  green  line.  The 
dotted  green  lines  represent  current  limits  of  the  motor.  The  current  results  are  saved 
in  Excel  file  idiq.xlsx. 
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4 


dq  Currents 


Figure  2.3;  Direct  and  quadrature  currents 
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4.  The  direct  and  quadrature  voltages  are  shown  here  in  Figure  2.4  and  are  similar  to  the 
currents.  The  direct  voltage  should  remain  around  0  and  the  quadrature  voltage 
should  contribute  to  the  requirements  of  the  selected  flight  profile. 


dq  Voltages 


Figure  2.4  (a):  Direct  and  quadrature  voltages 

The  Ud  voltage  does  not  vary  too  much  and  generally  stays  around  zero  as  seen  in  the 
zoomed  section  of  the  graph.  The  voltage  results  are  saved  in  Excel  file  uduq.xlsx. 
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dq  Voltages 


> 


T3 

C 

n: 

> 


60 


40 


20 


0 


?  -20 
> 


-40 


-60 


0.9 


0.95 


1.05 


1.1 


1.15 


Time,  t  [min] 

Figure  2.4  (b);  Zoomed  in  Section  of  Direct  and  Quadrature  Voltages 
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Power  Loss,  P  [kW] 

CL/  '■  ■* 


5.  The  copper  winding  inside  the  motor  contain  the  highest  source  of  heat.  The  transient 
heat  analysis  is  represented  here  with  peak  values  near  12  Watts  which  is  much 
higher  than  the  steady  state  power  loss.  The  results  are  saved  in  Excel  file  Pcu.xlsx. 

0.016 

0.014 

0.012 

0.01 

0.008 

0.006 

0.004 

0.002 

0 

0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 

Time,  t  [min] 

Figure  2.5;  Power  loss  in  the  motor  windings 


Power  Loss  in  Windings 
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6.  Of  great  interest  to  power  systems  engineers  is  the  power  in  (positive  values)  and  out 
(negative  values)  of  the  motor.  The  power  out  is  often  called  regenerative  power. 
Here  we  show  the  power  in  and  out  of  the  motor  in  Figure  2.6.  The  results  are  saved 
in  Excel  file  Pin.xlsx. 


Electrical  Power 

250  r 


-200  - 

-250 - ■ - 1 - ^ - ■ - 1 - 1 

0  50  100  150  200  250  300 


Time,  t  [min] 

Figure  2.6;  Electric  input  power  and  regenerative  power 
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Temperature,  r[C] 


7.  The  temperature  profiles  for  various  nodes  taken  from  different  areas  of  the  motor  are 
shown  below.  The  highest  temperatures  pertain  to  the  nodes  elosest  to  the  eopper 
windings.  The  ambient  temperature  was  held  at  22  °C.  This  flight  profile  showed 
only  a  O.d^C  temperature  inerease  over  5  minute  period.  The  results  are  saved  in 
temperature. xlsx. 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 

Time,  t  [min] 

Figure  1.1:  Temperature  profile  of  the  motor. 
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Part  III  -  Source  Code  (foc.m) 

%  Nonlinear  Dynamic  Modeling  and  Field  Oriented  Control 
%  of  Permanent  Magnet  (PM)  Motor 

g, 

o 

%  Written  in  SI  or  MKS  Unit  System. 

g, 

o 

%  Authors : 

%  David  Woodburn 
%  Dr.  Lei  Zhou 

%  Dr.  Thomas  X.  Wu 

clear  all;  close  all;  clc; 


%  Input  motor  parameters. 


motor Parameter 


g, _ 

o 

%  Load  EMA  mission  profile. 

g, _ _ 

o 

Mission  =  xlsread ( ' missionProf ile . xlsx ' ) ;  %  Load  into  mission  profile 
tN  =  Mission (:, 1 ).' ;  %  time  sequence 

xSN  =  Mission (:, 2 ).' ;  %  stroke  sequence 
FLN  =  Mission (:, 3 ).' ;  %  load  force  sequence 

%  xSN  =  xSN  *  0.0254;  %  If  stroke  is  in  inch,  needs  to  convert  to  meter. 

%  FLN  =  FLN  *  4.4482;  %  If  load  force  is  in  Ibf,  needs  to  covert  to  Newton. 

theta  meSN  =  xSN*Ncr*p/2  +  theta  melnitial; 
tauLN  =  FLN/ (Ncr*ef fDrive) ;  %  Load  torque  [N  m] 


g, _ _ 

o 

%  Set  up  simulation  parameters. 

5-  ___  ___  ___  ______________ 

o 

spc  =  10;  %  Samples  per  cycle  of  highest  frequency 

tauTh  =  0.01;  %  Thermal  time  constant  [s] 

%  Set  PI  coefficients. 

kp  =  0.1;  %  [] 

ki  =  0.0001;  %  [A  s/rad] 

%  Set  differential  time-stepping  weights. 
alpha_dt  =  0.5; 
beta_dt  =  1  -  alpha_dt; 

%  Set  motor  simulation  convergence  parameters, 
epsilon  =  0.0001;  %  Convergence  limit  [] 

Ncvg  =  50;  %  Maximum  number  of  convergence  iterations 


T  =  tN(end) 


tN  ( 1 ) ; 


T  could  have  been  redefined. 
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Nt  =  length (tN);  %  Initial  length  of  profile  arrays 

nfMax  =1;  %  Order  of  highest  frequency  harmonic  of  interest 

NtR  =  round ( spc* (nfMax*omega  meLim/ (2*pi) ) *T) ;  %  High-res,  time  steps 
dtEO  =  T/ (NtR-1);  %  Mean  hi-res  time  step 

omegaSMin  =  2*pi/0.002;  %  Min  sampling  frequency  [rad  sample/s] 
omegaSMax  =  2*pi/dtE0;  %  Max  sampling  frequency  [rad  sample/s] 

ntR  =2;  %  Hi-res  index  for  recording  [] 

NtETotal  =1;  %  Counter  for  the  total  number  of  simulation  points 
tTh  =  0;  %  Last  time  for  thermal  calculation  [s] 

%  Pad  arrays  (for  interpolation  referencing) .  Lengths  are  Nt+3. 
dtBar  =  mean (dif f (tN) ) ;  %  Note  that  dt  might  not  be  constant. 
tN  =  [tN ( 1 ) -dtBar,  tN,  tN (end) IdtBar,  tN (end) +2 *dtBar ] ; 

theta  meSN  =  [theta  meSN(l),  theta  meSN,  theta  meSN (end) ,  theta  meSN (end) ] ; 
tauLN  =  [tauLN(l),  tauLN,  tauLN (end) ,  tauLN (end) ] ; 

%  Initialize  states. 


Rs  =  Rs  room; 

q, 

o 

Initial  phase  winding 

resistance  [ohm] 

id  =  0 ; 

q, 

o 

Phase  d  current  [A] 

iq  =  0; 

q, 

o 

Phase  q  current  [A] 

did  dt  =  0; 

q, 

o 

Rate  of  change  of  id 

[A/s] 

diq  dt  =  0; 

q, 

o 

Rate  of  change  of  iq 

[A/s] 

LdO  =  LdSN (1) ; 

q, 

o 

Ld  at  zero  current 

LqO  =  LqSN (1) ; 

q, 

o 

Lq  at  zero  current 

Ld  =  LdO ; 

q, 

o 

Initial  Ld 

Lq  =  LqO; 

q, 

o 

Initial  Lq 

ud  =  0 ; 

%  Phase  d  voltage 

[V] 

uq  =  0 ; 

%  Phase  q  voltage 

[V] 

theta  me  =  theta 

meSN(l);  %  Start  up  (Mechanical  angle) *p/2  [rad] 

omega  me  =  0; 

%  (Mechanical  angular  frequency) *p/2  [rad/s] 

alpha  me  =  0; 

%  (Mechanical  angular  acceleration)  *p/2  [rad/s''2 

Qs  =  0; 

%  Heat  source  flux  (from  windings) 

Q  =  zeros (Nth, 1 ) ; 

%  Heat  flux 

NConverge  =  0; 

%  Counter  for  number  of  convergence  iterations 

%  Saturation  curve  points. 

rSat  =  0.10;  %  Saturation  ratio  [] 

il  =  iLim* ( 1-rSat) ;  %  Beginning  of  current  saturation  [A] 

i3  =  iLim* ( 1+rSat) ;  %  Ending  of  current  saturation  [A] 

ul  =  uLim* ( 1-rSat) ;  %  Beginning  of  voltage  saturation  [V] 

u3  =  uLim* ( 1+rSat) ;  %  Ending  of  voltage  saturation  [V] 

th2L  =  xMinLim*Ncr*p/2  +  theta  melnitial;  %  theta_me  max  limit 

th2R  =  xMaxLim*Ncr*p/2  +  theta  melnitial;  %  theta_me  min  limit 

thO  =  (th2R  +  th2L)/2;  %  Midpoint  between  theta  me  limits 

thlR  =  (th2R  -  thO) * (1-rSat)  +  thO; 

th3R  =  (th2R  -  thO) * (1+rSat)  +  thO; 

thlL  =  - (thO  -  th2L) * ( 1-rSat)  +  thO; 

th3L  =  - (thO  -  th2L) * ( 1+rSat)  +  thO; 

%  Set  threshold  to  show  limits  in  plots. 
limThreshold  =  0.8;  %  [] 
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%  Initialize  recording  arrays. 


%  Initialize  hi-res  arrays 
NtRO  =  round (NtR* 1 . 05 )  ; 
tR  =  zeros ( 1 , NtRO ) ; 
uaR  =  zeros ( 1 , NtRO ) ; 
ubR  =  zeros ( 1 , NtRO ) ; 
tauLR  =  zeros ( 1 , NtRO ) ; 
theta  meR  =  zeros ( 1 , NtRO ) ; 
omega  meR  =  zeros ( 1 , NtRO ) ; 


%  Add  5% 
%  Hi-res 
%  Hi-res 
%  Hi-res 
%  Hi-res 
%  Hi-res 
%  Hi-res 


buffer  [ ] . 

time  record  [s] 

phase  A  voltage  record  [V] 

phase  B  voltage  record  [V] 

load  torque  record  [N  m] 

theta_me  record  [rad] 

omega_me  record  [rad/s] 


%  Actual  theta_me  [rad] 

%  Actual  omega_me  [rad/s] 
%  Direct  current  [A] 

%  Quadrature  current  [A] 

%  Torque  angle  [rad] 


%  Initialize  low-res  arrays. 
tNR  =  zeros ( 1 , Nt) ; 
theta  meN  =  zeros ( 1 , Nt) ; 
theta  meN ( 1 )  =  theta  meSN(l); 
omega  meN  =  zeros ( 1 , Nt) ; 
idN  =  zeros ( 1 , Nt) ; 
iqN  =  zeros ( 1 , Nt) ; 
phiN  =  zeros ( 1 , Nt) ; 

LdN  =  zeros ( 1 , Nt) ;  LdN ( 1 )  =  LdO; 
LqN  =  zeros ( 1 , Nt) ;  LqN ( 1 )  =  LqO; 

udN  =  zeros ( 1 , Nt) ; 
uqN  =  zeros ( 1 , Nt) ; 

RsN  =  zeros ( 1 , Nt) ;  RsN(l)  =  Rs; 
tauMN  =  zeros ( 1 , Nt) ; 

PinN  =  zeros ( 1 , Nt) ; 

PcuN  =  zeros ( 1 , Nt) ; 

PfricN  =  zeros ( 1 , Nt) ; 

TthN  =  zeros (Nth, Nt) ;  TthN(:,l) 


%  Direct  inductance  [H] 

%  Quadrature  inductance  [H] 

%  Direct  voltage  [V] 

%  Quadrature  voltage  [V] 

%  Resistance  [ohm] 

%  Machine  torque  [N  m] 

%  Power  in  [W] 

%  Copper  loss  [W] 

%  Friction  loss  [W] 

=  Tth;  %  Temperature  [C] 


%  Run  through  time  steps. 


%  A  low  resolution  (low-res)  profile  is  analyzed  one  segment  (step)  at  a 
%  time,  breaking  each  segment  into  many  sub-segements  that  have  a  high 
%  resolution.  Interpolation  is  used  to  accomplish  this,  and  the  number 
%  and  size  of  the  time  steps  are  based  on  an  average  hi-res  time  step 
%  size  based  on  the  maximum  electrical  speed  of  the  motor.  During 
%  playback,  there  is  only  one  low-res  segment  (Nt  =  2)  and  all  the 
%  profile  points  are  hi-resolution .  At  that  point,  the  hi-res  time  step 
%  size  and  the  number  of  time  steps  are  based  on  the  pre-recorded  hi-res 
%  time  array.  Any  technique  used  for  interpolating  omega  on  the  low-res 
%  scale  will  be  yield  different  results  than  on  the  hi-res  scale  just 
%  because  of  the  scale,  even  if  the  same  technique  is  used.  However, 

%  using  a  perfect  N-point  derivative  (PNPD)  should  yield  fairly  accurate 
%  results  without  averaging. 

%  For  each  low-res  segment,  simulate. 
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for  nt  =  2 : Nt 

%  Calculate  number  and  size  of  steps  for  this  segment. 

dt  =  tN(nt+l)  -  tN(nt);  %  Update  low-res  segment  time-step  size. 

%  Calculate  variable,  hi-res  time  step  based  on  dthetaE/dt. 
dtheta  me  =  theta  meSN(nt+l)  -  theta  meSN(nt); 
dtheta_meAbs  =  abs (dtheta_me) ; 

omegaSAbs  =  spc*dtheta_meAbs/dt;  %  spc  =  samples  per  cycle 
omegaSAbs  =  omegaSAbs* (omegaSAbs>omegaSMin)  +  ... 

omegaSMin* (omegaSAbs<=omegaSMin) ; 
omegaSAbs  =  omegaSAbs* (omegaSAbs<omegaSMax)  +  ... 

omegaSMax* (omegaSAbs>=omegaSMax) ; 
dtE  =  2 *pi/ omegaSAbs; 

NtE  =  ceil(dt/dtE)  +1;  %  Update  number  of  sub,  hi-res  time  points. 
dtE  =  dt/ (NtE-1);  %  Update  hi-res  segments  time-step  size. 

%  Count  the  total  number  of  simulation  points. 

NtETotal  =  NtETotal  +  (NtE  -  1 )  ; 

%  Reset  micro-loop's  energy. 

PcuBar  =  0; 

%  This  is  used  to  find  the  average  power  in  a  micro  loop  which  is 
%  then  used  for  the  thermal  component. 

%  Simulate  this  low-res  segment  in  hi-res. 
for  ntE  =  2:NtE 


g, _ _ _ 

o 

%  Get  inputs  to  motor  simulation. 

q, _ 

o 

%  Variables  ending  in  "S"  represent  desired  values. 

%  Variables  ending  in  "P"  represent  next  step  values. 

%  Variables  without  these  are  this  step's  values. 

%  The  process  of  getting  the  inputs  to  the  motor  simulation 
%  transitions  from  this  time  step  to  the  next.  The  inputs  belong 
%  to  the  resolved  (convergent)  state  of  the  next  time  step. 

%  Therefore,  udP,  uqP,  tauLP,  thetaEP,  and  omegaEP  all  belong  to 
%  the  next  time  step.  The  motor  simulation  accepts  certain  of 
%  those  inputs  as  given  values  and  then  resolves  the  remaining 
%  states  of  the  motor  such  that  those  inputs  could  be  true. 

%  Determine  index  in  Nt  scale  from  index  in  NtE  scale, 
wt  =  ntE/NtE;  %  Progress  in  this  nt  segment  (0  to  1] 

%  Use  weighting  to  get  t,  thetaS,  and  omegaS . 

t  =  tN (nt) * ( 1-wt)  +  tN (nt+1 ) *wt;  %  [s] 

tauLP  =  tauLN (nt) * ( 1-wt)  +  tauLN (nt+1 ) *wt;  %  [N  m] 

theta  meS  =  theta  meSN (nt) * ( 1-wt)  +  theta  meSN (nt+1 ) *wt;  %  [rad] 

dt3  =  tN  (  (nt)  :  (nt+2)  )  -  tN  (  (nt-1)  :  (nt+1)  )  ;  %  [s] 

omega  me3  =  (theta  meSN (nt : (nt+2 ) )  -  ... 

theta_meSN ( (nt-1) : (nt+1) )). /dt3;  %  [rad/s] 

omega_meS  =  (wt  <=  0 . 5) . * (omega_me3 (1) . * (0 . 5  -  wt)  +  ... 
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wt)  +  . . . 
[rad/s] 


omega_me3 (2). *(0.5  +  wt))  +  ... 

(wt  >  0 . 5)  .  *  (omega_me3  (2)  .  *  (1 . 5  - 
omega_me3 ( 3 ) . * (wt  -  0.5));  % 

%  Get  theta  and  omega  errors. 
eTheta  me  =  (theta  meS  -  theta  me);  %  Mechanical  [rad] 
eOmega  me  =  (omega  meS  -  omega  me);  %  Mechanical  [rad/s] 

%  Get  iqS. 

tauM  iq  =  3*p/4* (lambdaPM  +  id* (Ld-Lq) ) ;  %  [N  m/A] 
iqS  =  iq  +  2/(p*tauM  iq) * ( I  * ( kp*eOmega  me/dtE  -  alpha  me)  +  ... 
c* (kp*eOmega_me) )  +  ki*eTheta_me/dtE;  %  [A] 

%  Set  ids. 

ids  =  0;  %  id  should  almost  always  be  zero  [A] 

%  Rate  limit  iqS  (reduces  transients) . 

diqdtS  =  (iqS-iq) /dtE;  %  Predicted  current  rate  [A/s] 

diqdtS  =  sign (diqdtS ) *min ( [didtLim, abs (diqdtS ) ] ) ;  %  Cap  diqdt  [A/s] 

iqS  =  iq  +  diqdtS*dtE;  %  Rebuilt  iqS  [A] 

%  Saturation  limit  abs  of  iqS. 

iqS  =  sign (iqS) *( (abs (iqS) <i3) * (abs (iqS)  -  (abs (iqS) >=il) * . . . 

(  (abs  (iqS)  -iLim*  (1-rSat)  )  ''2)  /  (4*iLim*rSat)  )  +  ... 

(abs (iqS) >=i3) *iLim)  ;  %  [A] 

omega  meP  =  omega  me  +  (kp*eOmega  me); 

%  Calculate  inductances  based  on  desired  currents. 

mLq  =  (NL-1) * (abs (iqS) -iLMin) / (iLMax-iLMin)  +  1;%  Unbounded  real  [] 
wLq  =  mod(mLq,l);  %  Real  number  [0:1)  [] 

nLq  =  mLq  -  wLq;  %  Unbounded  integer  [] 

LdS  =  (1-wLq) *LdSN (nLq)  +  wLq*LdSN (nLq+1 ) ;  %  [H] 

LqS  =  ( 1-wLq) *LqSN (nLq)  +  wLq*LqSN (nLq+1 ) ;  %  [H] 

%  The  above  is  equivalent  to  using  interpl,  but  much  faster. 

%  LdS  =  interpl ( iL, LdSN, abs ( iqS )) ; 

%  LqS  =  interpl ( iL, LqSN, abs ( iqS )) ; 

%  Calculate  new  voltages. 

udS  =  Rs*idS  +  LdS* ( idS-id) /dtE  -  omega  meP*LqS*iqS;  %  [V] 

uqS  =  Rs*iqS  +  LqS * ( iqS-iq) /dtE  +  omega  meP*LdS*idS  +  ... 

omega  meP* lambdaPM;  %  [V] 

%  Saturation  limit  abs  of  voltages. 

ud  =  sign (udS) *( (abs (udS) <u3) * (abs (udS)  -  (abs (udS) >=ul) * . .  . 

(  (abs  (udS)  -uLim*  (1-rSat)  )  ''2)  /  (4*uLim*rSat)  )  +  ... 

(abs  (udS)  >=u3)  *uLim)  ;  %  [V] 

uq  =  sign (uqS) *( (abs (uqS) <u3) * (abs (uqS)  -  (abs (uqS) >=ul) * . . . 

(  (abs  (uqS)  -uLim*  (1-rSat)  )  ''2)  /  (4*uLim*rSat)  )  +  ... 

(abs (uqS) >=u3) *uLim)  ;  %  [V] 


%  Seek  convergence  for  one  hi-res  time  step. 
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%  Assume  next  values  from  present  values.  Simulation  should  not 
%  assume  any  knowledge  of  "S"  values. 
idP  =  id  +  did_dt*dtE;  %  [A] 
iqP  =  iq  +  diq_dt*dtE;  %  [A] 

alpha  meP  =  alpha  me;  %  Angular  acceleration  [rad/s''2] 
omega  meP  =  omega  me  +  alpha  me*dtE; 

%  Run  dynamical  equations  until  convergence, 
for  nConverge  =  IrNcvg 


%  Get  next  inductances  by  look-up  table. 

mLq  =  (NL-1) * (abs (iqP) -iLMin) / (iLMax-lLMin)  +  l;%Unbounded  real 
wLq  =  mod (mLq, 1);  %  Real  number  [0:1)  [] 
nLq  =  mLq  -  wLq;  %  Unbounded  integer  [] 

LdP  =  (1-wLq) *LdSN (nLq)  +  wLq*LdSN (nLq+1 ) ;  %  [H] 

LqP  =  (1-wLq) *LqSN (nLq)  +  wLq*LqSN (nLq+1 ) ;  %  [H] 

%  The  above  is  equivalent  to  using  interpl,  but  much  faster. 

%  LdP  =  interpl ( iL, LdSN, abs ( iqP) ) ; 

%  LqP  =  interpl ( iL, LqSN, abs ( iqP) ) ; 

%  Get  next  magnetic  and  air  torques.  p/2  is  used  to  get 
%  cummulative  torque,  not  to  convert  to  electrical  units. 
tauMP  =  3*p/4*iqP* (lambdaPM  +  idP* (LdP-LqP) ) ;  %  [N  m] 

tauFricP  =  c*omega  meP*2/p;  %  [N  m] 

%  Get  next  angular  acceleration. 

alpha  meP  =  p/2/I*(tauMP  -  tauLP  -  tauFricP)  ;  %  [rad/s''2] 

%  Comes  from  tauM  =  tauL  +  I*alpha_m,  alpha  me  =  alpha_m*p/2. 


%  Get  next  angular  speed  and  angle, 
omega  meP  =  omega  me  +  ... 

(alpha_dt*alpha_me  +  beta_dt*alpha_meP) *dtE;  %  [rad/s] 
theta  meP  =  theta  me  +  ... 

(alpha  dt*omega  me  +  beta  dt*omega  meP) *dtE;  %  [rad] 

%  theta  meP  is  updated  here  because  some  parameters  might 
%  depend  on  the  value  of  theta  me. 


%  Get  next  currents.  Current  should  be  updated  last  since  it 
%  is  used  for  the  convergence  test. 

didP  dt  =  1/LdP* (-Rs*idP  +  omega  meP*LqP*iqP  +  ud) ;  %  [A/s] 
diqP  dt  =  1/LqP* (-Rs*iqP  -  omega  meP*LdP*idP  +  uq  -  ... 

omega  meP* lambdaPM) ;  %  [A/s] 

idP  =  id  +  (alpha_dt*did_dt  +  beta_dt*didP_dt) *dtE;  %  [A] 
iqPOld  =  iqP; 

iqP  =  iq  +  (alpha_dt*diq_dt  +  beta_dt*diqP_dt) *dtE;  %  [A] 


%  Check  for  electrical  convergence, 
if  abs(iqP  -  iqP01d)/iLim  <  epsilon  %  [] 
break; 

end 

end  %  end  nConverge 
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%  Add  convergence  iterations. 

NConverge  =  NConverge  +  nConverge; 

%  Update  heat  generated. 

Pcu  =  3/2*(iq^2  +  id^2)*Rs;  %  Next-step  source  heat  [W] 
PcuBar  =  PcuBar  +  Pcu/NtE; 


%  Simulate  thermal. 


if  t  -  tTh  >  tauTh 


%  Get  next-step  winding  resistance. 
Rs  =  Rs  room; 

dtTh  =  t  -  tTh; 
tTh  =  t; 


%  Calculate  temperatures. 

IC  =  Pcu/4;  %  winding  loss  (only  quarter  portion  is  used) 


%  Calculate  next-step  temperatures. 

TthP(l)  =  Tth(l)  +  (IC  +  (Tth (3) -Tth (1) ) /Rth (1) ) /Cth (1) *dtTh; 
TthP (2)  =  Tth (2) ; 

TthP(3)  =  Tth(3)  +  (IS  +  (Tth (4) -Tth (3) ) /Rth (2)  +  ... 

(Tth (10) -Tth (3) ) /Rth (12)  +  (Tth (1)  -Tth (3) ) /Rth (1) ) /Cth (3) *dtTh; 
TthP(4)  =  Tth(4)  +  ( (Tth (4) -Tth (14) ) /Rth (3)  +  ... 

(Tth (13) -Tth (4) ) /Rth (14)  +  (Tth (3) -Tth (4) ) /Rth (2)  +  ... 

(Tth (5) -Tth (4) ) /Rth (4) ) /Cth (4) *dtTh; 

TthP(5)  =  (Rth (5) *Tth (4)  +  Rth (4) *Tth (2) ) / (Rth (4) +Rth (5) ) ; 
TthP(6)  =  Tth(6)  +  ( (Tth (9) -Tth (6) ) /Rth (6)  +  ... 

(Tth (7) -Tth (6) ) /Rth (7)  +  (Tth (8) -Tth (6) ) /Rth (8) ) /Cth (6) *dtTh; 
TthP(7)  =  Tth(7)  +  (-IR7a  +  (Tth (6) -Tth (7) ) /Rth (7)  +  ... 

(Tth (12) -Tth (7) ) /Rth (9) ) /Cth (7) *dtTh; 

TthP(8)  =  Tth(8)  +  ( (Tth (6) -Tth (8) ) /Rth (8)  +  ... 

(Tth (13) -Tth (8) ) /Rth (10) ) /Cth (8) *dtTh; 

TthP(9)  =  Tth(9)  +  (IM  +  (Tth (6) -Tth (9) ) /Rth (6)  +  ... 

(Tth (10) -Tth (9) ) /Rth (11) ) /Cth (9) *dtTh; 

TthP(lO)  =  (Tth (3) *Rth (11)  +  Tth ( 9) *Rth ( 12 )  +  ... 

Rth (11) *Rth (12) *IW) / (Rth (11) +Rth  (12) ) ; 

TthP(12)  =  Tth(12)  +  (IBP  +  (Tth (14) -Tth (12) ) /Rth (13)  +  ... 

(Tth (7) -Tth (12) ) /Rth (9) ) /Cth ( 12 ) *dtTh; 

TthP(13)  =  Tth(13)  +  (IBB  +  (Tth (4) -Tth (13) ) /Rth (14)  +  ... 

(Tth (8) -Tth (13) ) /Rth (10) ) /Cth (13) *dtTh; 

TthP(14)  =  Tth(14)  +  (-IR14  +  (Tth (4) -Tth (14) ) /Rth (3)  +  ... 

(Tth (12) -Tth (14) ) /Rth (13) ) /Cth (14) *dtTh; 


Tth  =  TthP; 

end 


%  Update  states. 
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id  =  idP; 
iq  =  iqP; 
did_dt  =  didP_dt; 
diq_dt  =  diqP_dt; 

Ld  =  LdP; 

Lq  =  LqP; 

theta  me  =  theta  meP; 
omega  me  =  omega  meP; 
alpha  me  =  alpha  meP; 


g, _ 

o 

%  Record  hi-res  variables. 

q, _ _ _ _ 

o 

tHR(ntR)  =  t; 

udH(ntR)  =  ud; 

uqH(ntR)  =  uq; 

theta  meR(ntR)  =  theta  me; 

omega  meR(ntR)  =  omega  me; 

tauLR(ntR)  =  tauLP; 

ntR  =  ntR  +  1; 

end  %  end  hi-res  simulation 


%  Record  low-res  variables. 


%  Store  time. 
tNR(nt)  =  t; 

%  Store  motor  stroke. 

theta_meN (nt)  =  theta_me;  %  [rad] 

omega_meN (nt)  =  omega_me;  %  [rad/s] 

%  Store  currents. 
idN(nt)  =  id;  %  [A] 

iqN(nt)  =  iq;  %  [A] 

%  Store  inductances. 

LdN(nt)  =  Ld;  %  [H] 

LqN(nt)  =  Lq;  %  [H] 

%  Store  voltages. 
udN(nt)  =  ud;  %  [V] 

uqN(nt)  =  uq;  %  [V] 

%  Store  resistance. 

RsN(nt)  =  Rs;  %  [ohm] 

%  Store  motor  torque. 
tauMN(nt)  =  tauMP;  %  [N  m] 

%  Store  powers  [W] . 
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PinN(nt)  =  3/2*(iq*uq  +  id*ud) ;  %  Input  power 

PcuN(nt)  =  3/2*(iq"'2  +  id''2)*Rs;  %  Copper  loss  for  whole  motor 
PfricN(nt)  =  tauFricP* (omega  me*2/p); 

%  Store  temperatures. 

TthN(:,nt)  =  Tth; 

end  %  End  time  stepping. 

%  Unpad  profile  arrays. 
tN  =  tN (2 :  (Nt+1) )  ; 

theta  meSN  =  theta  meSN (2 : (Nt+1 ) ) ; 

tauLN~=  tauLN(2: (Nt+1) ) ; 

tR (ntR: end)  =  [ ] ; 

theta_meR (ntR: end)  =  [ ]  ; 

omega  meR(ntR:end)  =  [ ]  ; 

tauLR (ntR: end)  =  [  ] ; 

g,  _  _  _ 

o 

%  Show  statistics. 

g, _ 

o 


disp (' Dynamic  Modeling  and  Field  Oriented  Control  of  PM  Motor'); 


%  Calculate  stroke,  speed,  and  acceleration. 
xN  =  (theta  meN  -  theta  melnitial) *2/p/Ncr; 
vN  =  dif f (xN) . /dif f (tNR) ; 

NT  =  length (tNR) ; 

tNRmid  =  (tNR(l:NT-l)  +  tNR (2 : NT) ) /2 ; 
aN  =  dif f (vN) . /dif f (tNRmid) ; 


%  [m] 

%  [m/s] 

%  [s] 

%  [m/s^2] 


%  Calculate  and  Report 
PcuMean  =  mean(PcuN); 
tauMMax  =  max (abs (tauMN) ) ; 
tauLMax  =  max (abs (tauLN) ) ; 
FMax  =  tauLMax*Ncr*ef fDrive; 
max (abs (vN) ) ; 
max (abs (aN) ) ; 

=  min (diff (tHR)  ) ; 

=  max (diff (tHR) ) ; 


statistics . 

%  mean 


vMax  = 
aMax  = 
dtEMin 
dtEMax 
spaces  =  ' ' ; 
disp ( [ spaces 
disp ( [ spaces 
disp ( [ spaces 
disp ( [ spaces 
disp ( [ spaces 
disp ( [ spaces 
disp ( [ spaces 


winding 

[N  m] 

[N  m] 

[N] 

[m/s] 

[m/s^2 ] 

[s] 

[s] 


power  loss  [W] 


'Mean  winding  power  loss  =  '  num2str (PcuMean)  '  W  ']) 
'Max  motor  torque  =  '  num2str (tauMMax)  '  N  m  ' ] ) ; 

'Max  actuator  force  =  '  num2str (FMax)  '  N  ' ] ) ; 

'Max  velocity  =  '  num2str (vMax)  '  m/ s  ' ] ) ; 

'Max  acceleration  =  '  num2str  (aMax)  '  m/s''2  '  ]  )  ; 

'Min  electrical  time  step  =  '  num2str (dtEMin)  '  s  ']) 
'Max  electrical  time  step  =  '  num2str (dtEMax)  '  s  ']) 


%  Show  graphs . 


Stroke 


42 


figure ( ' Name ' , ' Stroke ' ) ; 

plot(tN/60,  xSN/0.0254,  tNR/60,  xN/0.0254); 

xlswrite ( ' strokeActual . xlsx ' ,  [tNR; xN] . ' ) ; 

ylabel ([' Stroke,  \it{x}\rm  [in]']); 

xlabel ( ' Time,  \it{t}\rm  [min]'); 

title (' Mechanical  Stroke  Following'); 

legend ( ' \it { x_{ desired] } ' , ' \it { x_{ actual } } ' ) ; 

%  Torque 

figure ( ' Name ' , ' Load  Force ' ) 
plot (tN/60, FLN) ; 

ylabel ('Load  Force,  \it { F_{ load] } \rm  [N]  '  )  ; 
xlabel (' Time,  \it{t}\rm  [min]'); 
title ( ' Load  Force ' )  ; 

figure ( ' Name ' , ' Load  Torque ' ) ; 
plot (tN/60,  tauLN) ; 

ylabel ('Load  Torque,  \it[ \tau_L} \rm  [N  m] ' ) ; 
xlabel (' Time,  \it{t}\rm  [min]'); 
title ( ' Load  Torque ' ) ; 

figure ( ' Name ' , ' Magnetic  Torque ' ) ; 

plot (tNR/60,  tauMN,  'g'); 

xlswrite ( ' tauM. xlsx ' ,  [tNR; tauMN]  .  ' )  ; 

ylabel (' Magnetic  Torque,  \it [ \tau_M} \rm  [N  m] ' ) ; 

xlabel (' Time,  \it{t}\rm  [min]'); 

title (' Magnetic  Torque'); 

%  Current 

figure ( ' Name ' , ' Current ' ) ; 

plot (tNR/60, idN, tNR/60, iqN) ; 

xlswrite ( ' i dig. xlsx' ,  [tNR; idN; iqN]  .  ' )  ; 

ylabel (' Currents ,  \it{i_d}\rm  [A]  and  \it{i_q}\rm  [A]  ' ) ; 

xlabel (' Time,  \it{t}\rm  [min]'); 

title ( ' dq  Currents ' ) ; 

legend ( ' \it { i_d} ' , ' \it { i_q} ' ) ; 

%  Voltage 

figure ( ' Name ' , ' Voltage ' ) ; 

plot (tNR/60, udN,  tNR/60,uqN); 

xlswrite ( ' uduq . xlsx ' ,  [tNR; udN; uqN]  .  ' )  ; 

%  xlswrite (' uduq . xlsx '  ,  [tHR; udH; uqH]  .  ' ) ;  %  for  high  resolution,  very  slow 

ylabel (' Voltages ,  \it{u_d}\rm  [V]  and  \it{u_q}\rm  [V] ' ) ; 

xlabel (' Time,  \it{t}\rm  [min]'); 

title ( ' dq  Voltages ' ) ; 

legend ( ' \it{u_d} ' , ' \it{u_q} ' ) ; 

%  Copper  Loss 

figure ( ' Name ' , ' Copper  Loss ' ) ; 

plot (tNR/60, PcuN/1000) ; 

xlswrite ( ' Pcu . xlsx ' ,  [tNR; PcuN] . ' ) ; 

ylabel (' Power  Loss,  \it { P_{ cu } } \rm  [kW]'); 

xlabel (' Time,  \it{t}\rm  [min]'); 

title (' Power  Loss  in  Windings'); 


43 


%  Electric  Power 

figure ( ' Name ' , ' Electric  Power ' ) ; 
xlswrite ( ' Pin . xlsx ' ,  [tNR; PinN]  .  ' )  ; 

PinNN  =  [ ] ; 
tNN  =  [  ] ; 
nO  =  1; 
hold  on; 

for  n  =  1 : length (tNR) -1 

PinNN  =  [PinNN,  PinN(n)]; 
tNN  =  [tNN,  tNR(n) ] ; 
if  PinN (n) *PinN (n+1 )  <  0 
PinNN  =  [PinNN,  0]; 

tadd  =  tNR (n) + (tNR (n+1 ) -tNR (n) ) *abs ( PinN (n) ) 
(abs (PinN (n) ) +abs (PinN  (n+1) ) ) ; 
tNN  =  [tNN, tadd]; 
if  PinN(n)  >  0 

plot (tNN/ 60 , PinNN, ' b ' ) ; 

else 

plot  (tNN/ 60, PinNN,  ' r ' ) ; 

end 

PinNN  =  0; 
tNN  =  tadd; 

end 

end 

PinNN  =  [PinNN,  PinN (end) ] ; 
tNN  =  [tNN,  tNR(end)]; 
if  PinN (end)  >  0 

plot (tNN/ 60 , PinNN, ' b ' ) ; 

else 

plot (tNN/ 60, PinNN, ' r ' ) ; 

end 

ylabel (' Power ,  \it{P}\rm  [W] ' ) ; 
xlabel ( ' Time,  \it{t}\rm  [min]'); 
title (' Electrical  Power'); 
hold  off; 

%  Thermal 

figure ( ' Name ' , ' Thermal ' ) ; 

plot  (tNR/60,  TthN  ([1,2, 3, 4, 5, 6, 7, 8, 9, 10,  12, 13, 14],;)) 
xlswrite ( ' temperature . xlsx ' ,  [tNR; TthN; ] . ' ) ; 
ylabel (' Temperature,  \it{T}\rm  [C]  ' )  ; 
xlabel (' Time,  \it{t}\rm  [min]'); 
title ( ' Temperature  Distribution ' ) ; 

toe; 
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Part  IV  -  Motor  Parameter  Code  (motorParameter.m) 


%  Define  primary  motor 
Rs  room  =  0.775; 
p  =  10; 

I  =  113.2*10" (-6) ; 
c  =  0.00001; 
lambdaPM  =  0.0793; 


parameters. (This  example  is  for  a  Danaher's  motor.) 

%  Phase  resistance  at  room  temperature  [ohm] 
%  Number  of  poles  [] 

%  Inertia  moment  [kg  m"2/rad] 

%  Friction  Coefficient  [kg  m"2/s  rad] 

%  PM  flux  amplitude  [Wb] 


%  Build  direct  inductance  Ld  and  quadrature  inductance  Lq. 
%  Load  current  and  inductance  arrays . 
iL  =  0:10:300; 


LdSN  = 


[4.46928937, 
055221239,  3. 
464592095,  3. 
096573444,  3. 
884038304,  2. 


4 
3 
3 
2 

2.739726553,  2. 
2.634980767,  2. 
2.553992315,  2. 
LqSN  =  [4.612083555 
3.518989894,  3. 
2.199031273,  1. 
1.568691069,  1. 
1.249202936,  1. 
1.054610808,  1. 
0.923114329,  0. 
0.828084863,  0. 


4.435673396 
895196891,  3 
34851825,  3. 
033981034,  2 
843532188,  2 
710652474,  2 
612739765,  2 
536488402,  2 
,  4.52697855 
09915783,  2. 
99252454,  1. 
470728535,  1 
192423128,  1 
017268661,  0 
896516968,  0 
808445183,  0 


4.334488 
.741845461 
250656076, 
. 978679457 
. 805950657 
. 683562881 
.592141217 
.519800427 
,  4.276907 
740369271, 
823051455, 
.386369068 
. 141458078 
. 98309393, 
. 871994934 
.789983231 


681,  4.201478319 
,  3.59804698, 
3.16810712,  . 

,  2.929048419, 

,  2.77189558, 

,  2.658379174, 

,  2.572607466, 

]  *  0.001;  %  [H] 
416,  3.919619295 
2.445203741, 

1 . 684390442, 

,  1.313170018, 

,  1.095896561, 
0.951930145, 

,  0.84926359, 

]  *  0.001;  %  [H] 


%  Define  thermal  parameters. 

Rth  =  [0.001389;  0.007697;  0.5;  0.001896;  16.0; 

9.1358;  16.97;  4.4581;  4.4;  5.933;  37.49;  2 
%  Thermal  Resistance  ["o  C/W] 

Cth  =  [38.85;  0;  82.3;  159.11;  0;  56.03;  14.99; 

6.72;  0.00176;  0;  13.18;  8.41;  62.56]; 

%  Thermal  Capacitance  [J  /  "o  C] 

Nth  =  length (Rth); 

Troom  =  22; 

Tth  =  ones ( length (Rth) , 1 ) *Troom; 

IS  =  0; 

IW  =  0; 

IBB  =  0; 

IBF  =  0; 

IM  =  0; 

IR7a  =  0; 

IRI4  =  0; 


0.08577; 
6;  2.6]; 

7.08;  .. 


Room  Temperature  ["o  C] 

stator  loss 
winding  loss 
rear  bearing  loss 
front  bearing  loss 
magnet  loss 


%  Get  current  min  and  max. 
NL  =  length (iL); 
iLMin  =  min(iL); 
iLMax  =  max(iL); 


%  Define  gear  train  coupling  ratio. 

Ncr  =  49.8728/0.0254;  %  Total  Coupling  Ratio  [rad/m] 
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strok;e=0 . 


effDrive  =  0.76;  %  Drive  train  efficiency  [] 

theta  melnitial  =  pi/6;%  Rotor  initial  (mechanical  angle*p/2 ) when 

kAir  =  0.00180297/ (2*pi)  ^1.927/ (2*pi*  (p*pi) '^0. 927); 


Define  simulation  limits. 


xMinLim  =  0; 
xMaxLim  =  4.05*0.0254; 
tauMLim  =  1.9256; 

PLim  =  688; 
uLim  =  165; 
iLim  =  19.2; 
didtLim  =  5*10''4; 
vLim  =  0.086; 


%  stroke  min  [m] 

%  stroke  max  [m] 

%  Motor  torque  limit  [N  m] 

%  Rated  output  power  [W] 

Maximum  instantaneous  voltage  allowed  per  phase 
Maximum  instantaneous  current  allowed  per  phase 
%  Maximum  current  change  rate  allowed  [A/s] 
%  Maximum  speed  allowed  [m/s] 


[V] 

[A] 


omega  meLim  =  vLim*Ncr*p/2 ;  %  Maximum  mechanical  angular  speed  *p/2  [rad/s] 
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